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Abstract. We obtain an elegant and useful description of the dynamics of 
Szekeres dust models (in their full generality) by means of "quasi-local" scalar 
variables constructed by suitable integral distributions that can be interpreted 
as weighed proper volume averages of the local covariant scalars. In terms of 
these variables, the field equations and basic physical and geometric quantities are 
formally identical to their corresponding expressions in the spherically symmetric 
» ' I _ LTB dust models. Since we can map every Szekeres model to a unique LTB model, 

j^j^ rigorous results valid for the latter models can be readily generalized to a non- 

spherical Szekeres geometry. The new variables lead naturally to an initial value 
formulation in which all scalars are expressed as scaling laws in terms of their 
values at an arbitrary initial space slice. Those variables also yield a significant 
simplification of numerical work, since the fluid flow evolution equations become 
a system of autonomous ordinary differential equations subjected to algebraic 
constraints containing the information on the deviations from spherical symmetry. 
As an example of how this formalism can be applied, we show that spherical 
symmetry is stable against small dipole-like perturbations. This new approach to 
the dynamics of the Szekeres solutions has an enormous potential for dealing with 
a wide variety of theoretical issues and for constructing non-spherical models of 
cosmological inhomogeneities to fit observational data. 



1. Introduction. 

The theory of General Relativity has been successfully employed over the years to 
describe self-gravitating systems, from the astrophysical up to the cosmological scale. 
Different phenomena need to be modeled by different solutions of Einstein's field 
equations. Assuming an idealized description, many self-gravitating systems can be 
examined by the well known Schwarzschild and Kerr solutions at the astrophysical 
scale, while the Friedmann-Lemaitre-Robertson-Walker (FLRW) models provide an 
adequate idealized framework at the cosmological scale. The extensive use of these 
simple exact solutions follows from the fact that the applicability of most of the 
thousands of known exact solutions to model astrophysical or cosmological systems, 
even as first order idealized approximations, is either impossible, or hard to justify 
and/or involves accepting unphysical constraints due to the symmetries characterizing 
these solutions [1]. 
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A lesser degree of idealization is provided by the well known and often used class 
of spherically symmetric exact solutions gcncrically known as the Lcmaitre-Tolman- 
Bondi (LTB) models [2, 3, 4]. These models allow us to examine non-linear effects 
in self-gravitating systems that have spherical symmetry or can be approximated by 
spherical configurations by means of exact analytic expressions or, at least, by using 
numeric but tractable methods. These models have been widely used recently to 
study the effects of cosmological inhomogeneities within the effort to explain cosmic 
observations without assuming the existence of a dark energy source (for a review see 
[5, 6]). Since spherical symmetry can be a strong and limiting constraint, even as 
a first approximation, it is important to look for exact solutions that consider non- 
spherical generalizations of LTB models, which provide a less idealized description but, 
at the same time, still allow for either an analytic treatment or a tractable numerical 
approach. The most general known class of solutions that generalize LTB models in 
this form are the Szekeres models, which (in general) admit no symmetries (no Killing 
vectors) and can be reduced to either LTB or FLRW models or the Schwarzschild 
solutions in the appropriate limits. 

The Szekeres solution was found in 1975 by Szekeres [7] and its first applications 
was the study of non-spherical collapse [8]. While there is a number of theoretical 
studies based on Szekeres models [9, 10, 11, 12, 13, 14, 15] (see reviews in [16, 5, 6]), 
only recently these models have gained more interest within the cosmological 
community as models to study light propagation [17, 18, 19, 20, 21], or structure 
formation [22, 23, 24]. Several new phenomena has been observed within the Szekeres 
model - for example the structures can evolve much faster than in the perturbed 
FLRW or Lemaitre-Tolman model [22, 23] or that two rays sent from the same source 
at different times to the same observer pass through different sequences of intermediate 
matter particles (as a consequence we should observe a drift of objects positions in 
the sky [25]). The Szekeres model was also used to provide a justification of existence 
of giant voids [26] - inhomogeneous alternatives to standard cosmological model. 

There arc two classes of the Szekeres model (sec [16] for detail): /?' 7^ and 13' ~ 
(using the Szekeres notation) or class I and II (using Goode-Wainwright notation [11]). 
Class 11/13' = solutions have not received much attention (though see [23, 24] for 
recent literature). Class 1/(3' 7^ solutions are better known because they are less 
idealized and easier to apply as models of inhomogeneities. This class of solutions has 
3 subclasses: "quasi-spherical" [12, 13], "quasi-planar" [14] and "quasi-hyperbolic" 
[15] (see [16, 5, 6] for a comprehensive review). These subclasses are not distinct as 
within one model there can be regions of quasi-spherical geometry followed by regions 
of quasi-hyperbolic geometry [15]. 

Since Szekeres dust models | are characterized by vanishing vorticity, 4- 
acceleration and magnetic Weyl tensor, they belong to the class of so-called "silent 
universes" , which are appropriate approximations to a general spacetime in the long 
wave limit [30, 31, 32]. Silent universes (in general) suffer from various problematic 
features, such as: linear instability and non-integrability of their spatial constraints 
in a fluid flow formalism [33, 34]. However, such problems do not arise for the case of 
Szekeres models and their higher symmetry subcases. As commented in [33, 34], there 
is a strong conjecture that Szekeres models may be the only spatially inhomogeneous 
self-consistent silent universes. 

I Szekeres models with uniform pressure where first examined by Szafron [27], see also [28]. For the 
use of q-scalars with nonzero pressure in the spherically symmetric case see [29] . 
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We aim in this paper to present a novel approach to the study of the Szekeres 
models of class I (/?' 7^ 0) that is based on new coordinate independent scalar 
variables ( "quasi-local" variables or "q-scalars" ) . These variables have been previously 
introduced in the study of spherically symmetric LTB models [35, 36, 37], looking at 
important theoretical issues, such as exact non-linear perturbations [35], averaging 
inhomogeneities [36, 38, 39, 40, 41], radial asymptotics [42], evolution of radial profiles 
[43], as well as a dynamical systems analysis [44, 45]. By extending the use of these 
variables to class I Szekeres models we can generalize these studies to non-spherical 
geometries. 

Since the q-scalars are defined in terms of suitable integral distributions of the 
local covariant fluid flow scalars (density, Hubble flow, spatial curvature), they become 
weighed averages of these local covariant scalars if defined as functionals (instead 
of functions). By expressing the local scalars as fluctuations of the q-scalars, these 
fluctuations together with the q-scalars and algebraic constraints linking them provide 
a scalar representation that completely determines that dynamics of the models, either 
in terms of analytic scaling laws or by means of autonomous evolution equations 
suitable for numeric work. 

It is a well known fact that Einstein's fleld equations reduce for all silent 
universes to a system of six fluid flow (or "1+3") evolution equations that contain 
no spatial derivatives, and thus can be treated (formally) as a system of six ODE's 
[46, 47, 48, 49, 50, 51]. If we use the standard fluid flow variables (as in the sources 
cited above) the spacelike constraints become a set of complicated non-linear PDE's 
on the spacelike coordinates. Since these PDE's must be satisfied al all times, they 
necessarily constrain initial data and make it harder to conduct the numeric work 
of integrating the evolution equations. However, if we use the q-scalars to construct 
a set of evolution fluid flow equations, not only we can handle these equations as a 
system of ODE's, but the PDE's that provide the spacelike constraints in the standard 
variables reduce to algebraic constraints, which are formally identical to those of 
LTB models, with the deviation from spherical symmetry entering through initial 
conditions. Hence, the new variables lead to a simplified dynamics for the Szekeres 
models and also allow for a better understanding of effects of their deviation from 
spherical symmetry. Evidently, dealing with algebraic constraints that are siipplied 
by means of initial data represents a valuable advantage over the traditional fluid flow 
variables used in [46, 47, 48, 49, 50, 51]. Also, the q-scalars may provide important 
information on key theoretical aspects, as well as illuminate the connection between 
the Szekeres and LTB models and with linear perturbations on a FLRW background. 

The section by section plan of the paper is as follows. We introduce in section 2 
the Szekeres models given in a parametrization that expresses their main quantities 
as formally identical to the corresponding quantities in LTB models. Field equations, 
covariant fluid flow scalars and their evolution equations are given in terms of this 
parametrization. A new set of quasi-local ("q-scalar") variables is introduced in 
section 3, together with fluctuations defined by comparing the local fiuid fiow variables 
with their associated q-scalars. We discuss briefly in section 4 the relation between 
the q scalars and averages, as well as with the decomposition of Szekeres scalars into 
monopole and dipole-like terms. Fluid flow evolution equations for the q-scalars and 
their fluctuations are obtained in various representations in section 5, showing how 
the deviation from spherical/planar/hyperbolic symmetry can be fully accounted for 
as initial conditions are specified. Initial conditions for these evolution equations are 
discussed in section 6. We introduce an initial value formulation in section 7, leading 
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to analytic solutions of the evolution equations in terms of the scale factors and scaling 
laws for the q-scalars and their fluctuations. In section 8 wc examine the regularity 
conditions to avoid shell crossings, which in the framework of this formulation can 
be stated as restrictions on the initial value functions. We show in section 9 how we 
can associate to each LTB model a 3 parameter class of Szekeres models by a simple 
transofrmation in the space of initial conditions. Each Szekeres model in this class 
is characterized by the same q-scalars as the LTB model, and its fluctuations can 
be obtained from the LTB fluctuations by a simple algebraic relation. In section 10 
we apply the initial value formalism to examine the stability of the deviation from 
spherical conditions by looking at the evolution of the dipole contribution in quasi- 
spherical Szekeres models close to spherical symmetry. In section 11 we summarize 
our results and provide a discussion of potential applications. We have included four 
appendices that complement the material covered in the article: scalar averaging 
in Szekeres models (Appendix A), covariant expressions for the q scalars and their 
fluctuations (Appendix B), the proof that functions of q-scalars are also q-scalars 
(Appendix C) and models with spherical or wormhole topologies (Appendix D). 



2. Szekeres models in the "LTB— like" parametrization. 



We begin with the spherically symmetric Lemaitre-Tolman-Bondi (LTB) metric in 
its standard representation: 



ds^ 



-dt^ + T^dr^ + i?2 (d^2 + gijj2 ^^^2^ 



(1) 



1-K 

where R = R{t, r), i?' = dR/dr and K = K{r). The corresponding field equations for 
a dust tensor and a comoving 4-velocity = Sq take the well known form 



2M' = 8TrpR^R' 



(2) 



where M = M(r) is the quasi-local or Misner - Sharp mass function and R = u'^VaR = 
dR/dt. It is straightforward to write the metric of the Szekeres model in a form that 
is similar to (1): 

ds2 = _dt2 + ^^dr^ + [dx2 + Ay^] , (3) 

where Y = Y{t, r, x, y) and £ = £{r, x, y) are given by 
R 



Y = 



P 



S 



x-Q 
S 



(4) 
(5) 



with S{r), P{r), Q{r) arbitrary functions, and e = 0,±1. The cases e = 1, 0, and — 1, 
respectively, correspond to quasi-spherical, quasi-planar, and quasi-hyperbolic models. 
The field equations take the same LTB-like form 



y2 = 



2M 



-K, 



where 



2M' = 87TpY^Y', 

M ~ K 



M = 



K 



(6) 
(7) 

(8) 
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Relevant covariant quantities take also the LTB-like form: the expansion scalar 
= VaU" and the Ricci scalar ^TZ of hypersurfaces t = const (orthogonal to u") 

2Y Y' 

e = i^ + ^. (9) 

(10) 

together with the shear and electric Weyl tensors 

~ Q 1 f Y' y\ 

0-a6 = V(aWb) - -^^o5 = S Sab, S = SabCT" ="3(^7"^)' (l^) 

= UcUdC'^'^'^ = W W = S„6 = _ ^ + |[ p, (12) 

where C"-"^'^ is the Weyl tensor, hab = 9ab + UaUb, Va = /i^Vb and Eab = Kb - ^aVb, 
with rja = y/hrrSa- The quantities (9)-(12) become identical to their LTB forms by 
replacing Y, K, M for R, K, M. 

The standard procedure to deal with the Szekeres model is to solve the 
Fricdmann like equation (6) by means of the following quadrature, which is equivalent 
to the LTB equation for R in (2): 

du 

Ju=o ^y2M/u - K 

where we eliminated Y^M^K in (6) in terms of R,M,K by means of (4) and (8). 
Above tkb{r) is another arbitrary function, the time locus of the big bang singularity, 
which adds to the available r-dcpcndent free functions S, Q, P, K, M. Considering 
that the metric (3) allows for an arbitrary rescaling r = r(f), so that any one of these 
functions can be eliminated by a suitable choice of the r coordinate, then the solution 
of (13) for a given choice of free functions determines a specific model in which the 
density follows from (7) and the remaining relevant quantities from (9)-(12). 

Since the scalars {p, 9, ^TZ, E, W} completely characterize the Szekeres models 
(considering (11) and (12)). an alternative to dealing with the dynamics of these 
models through the quadrature (13) is furnished by the "1+3" or fluid-flow approach 
[52], leading to the following scalar evolution equations: [52] 

p = -pQ, (14a) 
e = -|^-f /0-6E2, (146) 

E = -^GE- E^ + yv, (14c) 

>V = -eW-|pS + 3SW, (14d) 
together with the "Hamiltonian" constraint 

^ = ^-!^ + E^ (15) 
9 3 6 ^ ' 

and the spacelike constraints: 

Vb^a - lKQ,b = 0, Vb< - ^hy, = 0, (16) 
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Y'/Y 


3S ' 


Y'/Y 


as " 




— = 

3W ' 






Y'/Y 


Y'/Y 


3W 



which take the following form (component by component): 

(17a) 
(176) 

= 0, (17c) 

= 0, {17d) 

where and respectively denote d/dx and d/dy. 

Notice that the fluid flow system presented above is identical to that in the "1+3" 
fluid flow formalism [52] applied to LTB models [53, 54], save for the constraints 
(17c) and (17d) which contain the input from the "non-spherical" degrees of freedom 
through dependence on {x,y). Evidently, these constraints distinguish a Szekeres and 
an LTB model. Also, as a consequence of (17c) and (17(f), the quotient W/S only 
depends on {t,r), though W and S each depends also on {x,y). 

3. Quasi— local scalar variables. 

The hypersurfaces ^T[t] marked hy t = constant, whose Ricci scalar is given by (10), 
provide a natural time slicing for the Szekeres models. The ^T[<] in Szekeres models, 
in general, admit a foliation in terms of compact simply connected domains V C ^T[t] 
parametrized by the coordinates (r, x, y) (sec Appendix A §). For every scalar function 
A defined in V in an arbitrary '^T[t] we can generate the "quasi-local" scalar function 
(or "q-scalar") Aq : V ^ R dual to the local scalar A by means of the following 
correspondence rule: 

_ J^ATdV,,, ILjASY-Y'drdxdy 

' Iv^dV^P) JJJySY^Y'drdxdy ' ^ 

where = \/e — K and dV(p) is the proper volume element dV(p) = J^~^£Y^Y'drdxdy 
of the time slices ^T[t]. \\ 

Considering that (9) implies Q = [In(y^y)]', together with the forms of p and 
^TZ in (7) and (10), we obtain by applying (18) to ^ = p, Q, ^TZ the following closed 
analytic forms: 

Stt _ 2M _ 2M 

(20) 

!^ = A = A (21) 

§ The worldline r = cannot be contained in any domain V of quasi-hyperbolic models (e = —1). 
For quasi-planar and quasi-spherical models (e = 1) this worldline can be characterized as a special 
location where R = Y = T, = W = hold for all t. It is not a symmetry center in quasi— spherical 

models, since the latter are not spherically symmetric. 

I The term "quasi-local" follows from the integral definition of the quasi-local Misncr-Sharp mass- 
energy function in spherical symmetry [35]. The relation between (18) and the average of A with 
weight factor T over the domain V is discussed in section 4.1 (see also [36, 40, 41]). 
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so that the q-scalars pg, 6g, "^TZg only depend on (t,r), even if their "local" 
counterparts p, Q, ^TZ depend (in general) on all four coordinates. % The relation 
between (19)-(21) and the weighed averages {p)q, (6)g, {'^T^)q obtained from treating 
(18) as a functional is discussed in Appendix A (see section 4.1 and [36, 40, 41]). 

By differentiating (19), (20) and (21) with respect to r and applying (7), (9) and 
(21) we obtain 

p; = ^(p-p,), e; = ^(e-e,), ^n'^ = ^-^{^Ti-'n,), (22) 

which implies, by means of (9). (11) and (12), that the scalars associated with the 
shear and electric Weyl tensors are expressible as fluctuations of Q and p with respect 
to their q-scalar duals 9g and pg-. 

1 47r 

E = --(e-e,), w = -—ip-pg). (23) 

In order to connect local and q-scalars, we define the relative fluctuations of p, 6, '^TZ 
as 

Ag SY'IY' ^' ' ' ^ ^ 

where we used the derivation rule (22). Since pg, Og, ^TZg are independent of {x,y), 
we have from (24) for A = p, &, ^TZ: 

= -A(^) = A(^) (25a) 

^(.) ^ _^(.) 0^ ^ ^(.) (256) 

where we used Y' /Y = R' / R — £' /£ from (4). As shown in Appendix B, the 
fluctuations (24) are coordinate independent quantities related to curvature invariants. 



4. Some appealing properties of the q— scalcirs. 

4.-1- The q-scalars as averages. 

If we use the correspondence rule (18) to define functionals (not functions) then we 
can interpret the q-scalars as proper volume average distributions {A)g (see Appendix 
A) on a domain V with weight factor 

1/2 



(26) 



which in spherical symmetry is a scalar invariant that reduces to the "7" factor in 

the Special Relativity limit and to total (or "binding" ) energy in the Newtonian limit 
[55]. The difference between the q-scalar and the q-average is subtle: for any given 
domain T> the q-average associates the number {A^g-D to the whole domain D, while 
the q-scalars are pointwise local functions 7? R, hence both are equal only at the 
boundary dD of each domain (which is marked by r constant). As a consequence, 
both satisfy the same local derivation rules at the boundary of an arbitrary D, but 

^ While the q-scalars pq, Oq, ^TZq axe coordinate independent quantities (see Appendix B), the 
function M given in (7) is not an invariant scalar of the Szekeres models (as its equivalent M is for 
spherically symmetric spacetimes). 
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behave differently under integration over the whole V. An elaborate discussion on the 
relation between q scalars and q-avcragcs is given in [36, 40, 41] for the spherically 
symmetric case (though it applies to the Szekeres case). 

As shown in [56], the standard proper volume average (scale factor = 1) of a 
local scalar A in all Szekeres models is independent of {x,y). Wc show explicitly 
in Appendix A that this is also the case for the q-scalars. Unless explicitly stated 
otherwise, we will consider throughout the article the Aq obtained from (18) as new 
local functions whose correspondence rule is the same as that of the weighed average 
functional 



4-2. Decomposition in terms of a monopole and a dipole. 

It is important to remark that relevant physical and geometric scalar quantities in 
Szekeres models, such as density, expansion rate, etc, can be decomposed as a sum of 
a pure "radial" part and a "non-radial remainder" + , i.e. 

A = + A„„.„d, Aad = A\£,=o, (27) 

where A^^^ only depends on t and r. As shown in [56], in computing the standard 
proper volume average of A (which is (18) with T = I) for quasi spherical models 
(e = 1) the dipole cancels out, so that {A) = {A^^^}- We prove in Appendix A that 
this also holds for the q-average with weight factor T ^ 1 with the "angular" part 
canceling out also in the quasi-hyperbolic and quasi-plane cases, i.e. we have for every 
domain X'[r] bounded by a surface r = constant: 

, , , _ X- ^ £y^y'drdxdy _ J A^,R^R'dr _ ^ 

J J J SY^Y'drdxdy ~ ^ J R^R'dr " <^)«^H' ^^8) 

where we are assuming that these integrals are bounded for a given domain contained 
in the slices ^T[t\ in models that are not quasi-spherical. As a consequence, wc can 
think of q-scalars as providing at every domain X>[r] the q-average of the monopole 
"radial" term in the decomposition (27). 



5. Evolution equations for the quasi— local variables. 

Considering (23) and (24), the five scalars {p, ©, S, W, ^TZ] that characterize the 
fluid flow dynamics of the Szekeres models are expressible in terms of the scalar 
representation {Aq, A^^^} for ^ = p, 9, ^Tl: 

p = Pq{l + /\^py), n=nq{l + A^^'>), /C = /C,(l + A('C)), (29a) 

A-TT 

i: = -'HgA^^\ W = -—pgA(p\ (29&) 
where we are using (and will use henceforth) the notation 

e 

U^-, /C^— . (30) 
so that A(®) = A(«), A('^) = AC^). 

+ In quasi-spherical case this "remainder" has the structure of a dipole, while in the quasi-hyperbolic 
case it is a pseudo-spherical equivalent of a dipole. 
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5.1. Evolution equations and constraints. 

Eliminating the local scalars {p, "H, S, W} in terms of {pq, Hq, A^^^ A^^^} by means 
of (29a) (29&) and substituting into the 1+3 system (14a)-(14(i) yields the following 
system of autonomous evolution equations 

Pq = - 'iPqHq, (31a) 

47r 

Hq =-nl-—pq, (316) 

A^P^ = -3{l + A^f>^)HqA^^\ (31c) 
A(«) = - (1 + 3A('^)) ^gA(«) + 1^(A(«) - A(^)), (31d) 

while substitution of (29a)-(29&) in (15) and (16) yields the algebraic constraints 

'Hq = ^Pq-K^q, (32) 

2 A(«) = Qq + {1 - nq) ACC) , (33) 
where we have introduced the following q scalar analogous to a FLRW Omega factor 



nq^-^, ^q-l = ^. (34) 

whose corresponding fluctuation is (from (24)) 

^ ^(p) _ 2a(«) = (1 - 0,)(A('') - A('=)), (35) 

and its local dual is 17 = 17,(1 + A^^^). As we show in Appendix C, all scalars 
expressible as functions of q-scalars are also q-scalars and comply with (18) and (24). 

We remark that the constraints (17a)-(176), associated with the "radial" 
derivatives of S and W, reduce to the first two identities in (22), while the constraints 
(17c)-(17d), associated with the "non-radial" derivatives of S and W, reduce to the 
identities (25a)-(256), which are valid for all t. While the four constraints (17a)-(176) 
of the system (14a) (14rf) are differential equations involving spatial gradients, the 
constraints (32)-(33) of the system (31a)-(31d) are algebraic relations that will hold 
for all t once they hold in an initial slice t = to, and thus can be used to set up the 
initial conditions in terms of the coordinates (r, x, y). 

As a consequence, when integrating the system {31a)-{31d) the only spacelike 
constraints that need to be taken care of are (32)-(33), which are algebraic relations 
that can be used to set up the initial conditions in terms of the coordinates (r, .T,?y). 
The evolution equations (31a)-(31(f) can be treated then as a system of ordinary 
differential equations constrained by a 3-parameter set of initial conditions satisfying 
(32)-(33) at an initial time slice. 



5.2. Alternative scalar representations. 

The evolution equations (31a)-(31(i) are given in the representation {pq, Hq, A'-''\ A^^^}. 
However, the q-scalars Aq = {pq, JCq, Hq, fig} are interrelated by the two constraints 
(32) and (34), while their fluctuations A^^' = {A^''), A^^), A(«), A^")} are inter- 
related by the two constraints (33) and (35). Hence, it is suflicient to select any 
representation {Aq, A^^^} made by any two of the Aq and any two of the A*^^^ (with 
B ^ A in general) to determine completely the dynamics of LTB models through 



A novel approach to the dynamics of Szekeres dust models 



10 



evolution equations analogous (and equivalent) to (31a) (31 d). Considering the re- 
lation between Hg, Qg and cosmological observable parameters, a useful alternative 
representation given by {fig, 'Hg, A^^\ A^^^} follows by using (32)-(35) to eliminate 
pq and A^^^ in terms of fig and A^^^ , leading to 



Hg = -l^l + lngjul (36a) 

iig = -ngii-ng)nq, (366) 



A(«) = - 



(l + 3A(«)) A(«) + ifig (a(«) + A(")) 



Hg, (36c) 



17gA(") + (rig - 3A(") - l) aC")] Hg, (36d) 

As the system (31a)-(31d), these evolution equations can also be treated as a system of 
autonomous ODE's subjected to the same algebraic constraints. However, equations 
(36a)-(36d) are more suited to be used to generate an dynamical systems study for 
Szekeres models as has been done with the LTB model with and without a cosmological 
constant [44, 45]. 

6. Initial conditions. 

Since {pg, Hg, Kg, O^} do not depend on {x,y), they become r-dependent functions 
{PqOi Hgo, JCgo, ilgo} lu au initial slice ^T[to] for an arbitrary to (the subindex o will 
denote henceforth evaluation at t = to). On the other hand, the initial value forms for 
the relative fluctuations {Ao'\ Aq'^', Aq^\ ^o^^} will depend on {r,x,y) through the 
function £ in (5), whose gradient £'/£ enters in the definitions of these fluctuations 
in (24), which are valid for all t. Specifying £ requires prescribing the three arbitrary 
r-dependent functions S, P, Q. 

Because of the constraint (32), we can choose any two of the four functions 
{pgo, Hqo, ICgi, flgo} as initial value functions. These two initial value functions, 
together with S, P, Q, are sufficient to determine the whole set of initial conditions 
for evolution equations like (31a)-(31(f) or (36a)--(36(i) by means of (24), (32), (33) 
(34) and (35). It is useful to fix the r coordinate also in the initial slice by the choice * 

i?o=r, so that Yo = ^, ''^ = \-l^. (37) 

As an example, if we choose the set of initial value functions {pgo, tCqo, S, P, Q}, the 
remaining initial value functions needed to integrate (31a)-(31c?) are 



-172 Stt Sirpgo 

"5 '^qO 



A^^^ = A = p,IC, H, n, (386) 

with the 5q^^ being the LTB initial fluctuations obtained from (24) at t = to for LTB 
models {£' = 0, Y = R): 

4^^ = ^^, si-^ = '-^, (39a) 

6 PgO O l^qO 

* This coordinate choice is not appropriate if the slices ^T[t] have spherical (S^) or wormhole (S^ x M) 
topologies. In these cases, Ro{r) must have two zeroes or no zeroes. We look at these cases in 
Appendix D. 
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l-qO r(/C) 



(1 - - 5, 







(396) 
(39 c) 



3 HqO 
^ 3 flqO 

where we used (24), (32), (33), (34) and (35). For the alternative system (36a)- 
{36 d) the appropriate choice of initial value functions is furnished by the set 
{Tiqi, flqo, S, P, Q}. This is a practically useful alternative since these functions can 
be related to the Hubble an Omega parameters a,t t = to. In this case, the remaining 
functions are 



2 t~\ 
a{p) _ "0 + 

° 1 - r£'/£ 



qO 



nio i^qo - 1), 



° 3(1- r£'/£) ' 



[nqo/{nqo - 1)] <5^''^ + 25, 



(40 a) 



(406) 



1 - r£'/£ 

The relation between these initial functions and the standard free parameters M and 
K follows readily from (19), (21) and (37): 



(41) 



2M = —pgo = 0,0-^^0 r\ K = ICqo = {^Igo - m^o 

while the bang time can be obtained as a function of the initial value functions 
(for example pqQ, JCqo) from the solutions of (6) (sec section 6.2). 

In order to prescribe the "non-spherical" part of the initial conditions it is useful 
to transform the {x,y) coordinates of (3) by means of a stereographic projection of 
polar coordinates [12, 15], which takes the following form: 

{S'cot(6'/2)cos((/)), Scot (61/2) sin(^)} for e = 1, 



{x-P,y-Q}={ 



{S{e/2)cos{(j)), S{e/2)sm{(P)} 



for 



0, 



(42) 



{5coth {e/2) cos((/)), S'coth (9/2) sin((/))} for e = -1, 

and leads to a non-diagonal metric because the transformation also depends on r 
through S, P, Q (see [13, 22]). However, this has no consequence for our purposes 
which is to parametrize an appropriate domain of initial conditions by means of angular 
coordinates {9, (p) with finite ranges and clear geometric interpretation. Applying (42) 
to (5), we can rewrite £'/£ as [12, 15] 

- [S' cos 61 + sin 61 (P' cos (t> + Q' sin 0)] /S 



£' 

j = { - [S' 6 {P' cos (f> + Q' sin (f>)]/S 



for e 
for e 



0, (43) 



- [S' cosh 9 + sinh 9 (P' cos </> + Q' sin cj))] /S for e = - 1, 
so that the particular case with axial symmetry (e = 1 and Q, P constants) yields 

(44) 

which clearly has a dipolar formjj. 



£' S' , 
-=--cos^. 



tl Note that this is only one particular representation of the axially symmetric case. Other forms that 
include P and Q functions are also possible. This representation, however, has the simplest structure. 
For other parameterizations of the axially symmetric cases see [5] . 
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Since practically all relevant expressions associated with the Szekeres models are 
formally identical with the corresponding LTB expressions, save for the presence of 
the term 1 — r£' /£ (which is multiplied or added), then it is very useful to employ 
(in numerical and analytic computations) its form (43) in angular coordinates, or (44) 
for the case with axial symmetry, even if these quantities have been obtained in the 
standard {x, y) coordinates. However, we will use angular coordinates only to calculate 
£' /£ in setting up initial conditions, with the time evolution obtained from systems 
like (31o)-(31d) or (36a)-(36d). 

We emphasize that the initial conditions (38a)-(38&) (or (40a)-(40&)) have 
been obtained by applying the algebraic constraints (32)-(33) of the system, while 
the constraints involving spatial gradients of the fluid flow system ((17a)-(17(i)) 
correspond to mathematical identities satisfied by the fluctuations A^^^ ((24) and 
(25a)-(25&)) that are satisfied for all t, and thus are satisfied &t t = t^. As a 
consequence, given any set of initial value functions (which must include S, P, Q), 
the dynamics of any quasi-spherical Szekeres model can be fully determined by the 
system (31a)-(31d), which can be integrated numerically by means of techniques used 
for autonomous ODE's. 



7. An initial value framework. 



Together with the evolution equations described before, we can examine the dynamics 
of Szekeres models by analytic expressions given in terms of scaling laws with respect 
to a given set of initial conditions (as in [42, 43] with LTB models). Evidently, these 
scaling laws are the analytic solutions of the evolution equations. 
We define the dimensionless scale factors 



r = 



Yo Ro 
Y'/Y 



r - r£'/£ 



with: r = 



R'/R 



= 1 + 



ra 



(46) 



Y^/Yo l-r£'/£' ' R'JR^ 

where we used the coordinate choice (37). The metric (3) takes the FLRW-like form 



'dt' 



cl//) 



£2 



(47) 



7.1. Scaling laws. 

Since M = M{r) and K = K{r), considering (19), (21), (30) and (41), and comparing 
the expressions for p and /C in (29a) with the expressions for p and ^TZ in (7) and (10), 
together with (8) and (24), we obtain the following scaling laws in the representation 

{p, /c, Mp\M'^)} 



Pq 



PqO 
n3 ' 



/c, 



1 + a(^' 



+ ac^) = 



1 + A^ 
f 

2/3 + A, 



1 + S^f'^ - r£'/£ 
r - r£'/£ ' 



(IC) 



m{l-r£'/£) 
T-r£'/£ 



(48 a) 
(486) 

(48c) 
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where f and F are defined in (46) and S^^^ , Jq'^^ are the LTB initial fluctuations given 
by (39a)-(396). The scahng laws for the remaining q-scalars and fluctuations follow 
readily from (20), (32), (33), (34) and (35): 



nj-i 

TLft n 



qO 



n 



2mgo 



where 2mqo 



2mqa - JCqoa 
Tr/3)pqo and 



' 2mqo-JCqoa' 



a(«) 



qO 



+ (1 - 0,o)a 



(1-r) 



2[nqo + {l-nqo)a] (T-rS'/S) 



(1-U 



qO) 



1(1 -r) 



(49) 

(50) 

(51) 
(52) 



[nqo + {1 - nqo) a] {T~r£'/£) ' 

It follows from these scaling laws that for all models pq, fCq, T-Lq diverge as a — > 0, 
while Qq 1 holds in this limit. On the other hand, all the A^"*^ diverge as F ^ 
(or F — >■ r£' /£), which marks a shell crossing singularity if it occurs for a > 0. 

7.2. Analytic solutions. 

Implicit analytic solutions of the quadrature (13) in the representation 
{p, K, A'-p\ A^'^^} depend on the sign of ICqo, leading to parabolic (/Cgo = 0), hy- 
perbolic (/Cgo < 0) and elliptic (/C^o > 0) models 

2a3/2 



Parabolic, 
Hyperbolic, 

Elliptic, 



t-i>,H=-^^=, (53) 

(54) 

t-t^^=( (55) 



t ^bb — 



3^2mqo ' 
Zhjap a) 

Ze{ao a)//3o 



expanding phase : Tiq > 0, 



[27r — Ze{ao a)] //3o collapsing phase : Hq <0, 
where ao = \lCqo\/mqo, /3o = \K^qo\^^'^ /mqo and Zh and Zg are 

u ^ Zh{u) = u^l"^ (2 + ufl'^ - arccosh(l + u), (56a) 
u Ze(u) = arccos(l - u) - u^/^ (2 _ ufl"^ . (566) 

Since a = 1 at t = to and a = at t = f^b and t = t^^^y (for elliptic models), we obtain 
the bang and collapse time functions 
2 

tbh = to — - — 1^ Parabolic, (57a) 



to 



ibb = to — 



3y^2m5G 
Zhjap) 

Ze{ao) 



Hyperbolic, 



tcoll — tb 



2tt 
Tp 



tp + 



27r-Ze(Qo) 



Elliptic, 



(57&) 
(57c) 



Po PO po 

Notice that it is possible to use tbb as an initial value function. This requires providing 
a specific functional form for tbb and (say) pqp, then we can flnd ICqp by solving (576) 
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or (57c) numerically. Also, the case with a simultaneous big bang transforms (57&) 
and (57 c) into constraints so that given mgo = (47r/3)pgo we find /C^o by solving them 

with = ^bb^ = constant. 

Note that the above relations do not depend on (x, y) variables and are the same 
for the quasi-spherical, quasi-hyperbolic, and quasi-plane cases. The dependence on 

{x,y) variables enters only via initial conditions Aq'''' and/or Aq'^-'. 

By implicit derivation of the solutions (53), (54) and (55) we obtain the following 
form for T given by (46): 



P„ab„,ic r^, + ^i"--0^^^. (58.) 

Hyperbolic and elliptic 

f = 1 + 3(A('') - A^'^)) - 3 (a(^) - ^Af Hqit - i,,) - Y^^^' (58&) 
where in the hyperbolic and elliptic case t^h is given by (576) and (57c) and 

rC A(^>-Af) / (,) 3 (^)\ 
3{l^r£'/£) - —R;, ~ ) " ^^^^ 

follows from deriving (576) and (57c) and using (76). All local scalars {p, JC, %, S, W} 
can now be given as functions of a and the initial value functions by applying the 
scaling laws (48a)-(51) to the relations (29a)-((296), with f given by (58a) or (586). 

7.3. Scaling laws and analytic solutions in other representations. 

7.3.1. The representation {il,, "Hg, A(^\ A^'^^} . The scaling laws and analytic 
solutions follow from the expressions derived before by replacing and /C^o with fi^o 
and "Hgo by means of (40a) (406), so that parabolic, hyperbolic and elliptic models 
now correspond respectively to fi^o — 1 = 0, fi^o — 1 < and fi, — 1 > 0. The scaling 
laws (48 a), (49) and (50) take FLRW forms 

Stt f^qo'H^o ^ (^^gii - l)Hyo 

-irPq = 5-^, = 2 ^' (60a) 



a/2 _ a/2 

il-q — 'T-gO 



(606) 



while the fluctuations A^"^^ take the same form as in (486), (48c), (51) and (52). The 
analytic solutions (53)-(55) and the forms of tbb in (57a)-(57c), as well as previous 
expressions for T and t'^^^^ have the same forms with ao and /3o given by: 

7.5.^2. The representation {Vlq, Hq, A^^^} . If we keep a as scale factor 

then the scaling laws for the q-scalars and their fluctuations are the same as in 
the representation {Vlq, Hq, A^p\ A^'^^} above, with Aq'\ Aq^^ expressed in terms 
oiA^^\A^'^ by (406). 
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Another possibility follows by using Og as scale factor by eliminating a and Hq 
in terms of Hqo, flqo, ^q by means of 



flqO (1 - »g) 

Og (1 - 



qO) 



T-Lq — T-LqO — 



1 - n. 



qO 



3/2 



(62 a) 



(626) 



The analytic solutions of the quadrature (13) become trivial for the parabolic case 
(Og = 1), while for the hyperbolic and elliptic case they take the form 

W-Wo 



hyperbolic 
elliptic 

t-to-- 



< Ogo < 1 : 



t-to 



n 



qO 



figO > 1 

w- 



Wo]/K 



qO 



[7rOgo(f^gO - 1) 



(expanding phase), 

-3/2 _ ^ _ ^ 1 



qO 



(collapsing phase), 



where the functions W and Wo 
eof^go 2y'£o(l - 



W\ 



take the form 



W 



211 



l^g) 



0„ 



A 



Wq = 



|1 — ^qo\ 

where Eq = 1, 



1 - 



qO 



11/2 



A 



2 

^qO 



- 1 



2|l-ng0 

A — arccosh correspond to the hyperbolic case (Ogo < 



-1, ^: 



arccos to the elliptic case (figo > 1)- The bang time is given by 



tbh =to— 



'HqO 



(63) 

(64) 

(65) 

(66) 
1) and 

(67) 



while r and the gradient of the bang time follow from (586) and (59) with Hq given 
by (626) and Ag'''', Aq'^-' eliminated on terms of Aq^-*, Aq^'' by means of (406). 

Using the functions Hq, fig and their initial values as variables of the scaling laws 
can be very useful for future work in applying Szekeres models to fit observations, 
as Hq, fig provide an appealing generalization of FLRW cosmological observational 
parameters (they reduce to these parameters in the FLRW limit). In fact, LTB void 
models used to fit observations are often parametrized in terms of Hq, fi,, exactly 
defined as in (62a) and (606), which are introduced as ansatzes [57, 58, 59, 60, 61, 
62, 63, 64]. Since Hq is a QL variable and fig is the ratio of QL variables, and the 
latter are Szekeres variables independent of the "non-radial" coordinates x, y, then it 
is expected that they exactly correspond to (and satisfy the same relations) as their 
analogous LTB variables. 



8. Regulcirity Conditions. 

As with LTB models, Szekeres models admit two types of curvature singularities: a 

"central" singularity associated with a — and t — t^^i, or (in elliptic models) t — ^con , 
and a shell crossing singularity T{t, r, x, y) = 0. Notice that the "non-radial" variables 
{x, y) play no role in determining the locus of o = 0, but they are involved in the shell 
crossing singularity (and in the conditions to avoid it). Also, it is important to remark 
that a = Y = Q and F = F' = 0, but the converses of these implications 
are not true. 



A novel approach to the dynamics of Szekeres dust models 16 

The condition to avoid a shell crossing singularity is given by F > and by 
using (58a) and (586) it can be expressed as conditions on the initial value functions 
(the generalization of Hellaby-Lake conditions in LTB models [65, 66]). Particular 
conditions depend whether we have the parabolic, hyperbolic and elliptic cases: 

Parabolic models. It is evident from (58a) that the necessary and sufficient 
condition for F > is given by: 

-1 < A^^^ < 0, (68) 

where we used the relation between Aq''-' and Sjf ^ in (39 a) together with Sjf'' = 

HqQ rt^^^ . Notice that the condition involving Aq'''' is equivalent to M' > and 

Po = Pgo[l + ^o'^] > 0) where M is defined by (8) and po is the initial local 
densityff. Conditions (68) reduce to their equivalent forms for parabolic LTB 

models {£' = 0) [37, 42, 43, 65, 66]. 

Hyperbolic models. We look at the form of F in (46) in the following asymptotic 
limits along comoving worldlines: 

• a-^0, Uq^oQ, nq{t-'t^^)^l + 0{a) 

f « 1 + A^^) - rUq C/(l - t£'/£), (69) 

• a^oo, Hq^Q, 'Hq{t-t^^)^l + 0{^) 

f«l + ^A('^). (70) 
Hence, necessary and sufficient conditions for F > are given by: 

a(^)>-1, a('=)>-| _^<o, (71) 

As in the parabolic case, the condition given in terms of /S.\!^^ implies that M' 
and Po ssc non-negative, while the condition given in terms of Aq'^' implies that 
K' >Q (notice that ^ < and /Cgo < hold for hyperbolic models). Conditions 
(71) reduce to their forms for LTB models if f = so that Aq^-* = d^'^ and 
Aq'^' = 4'^' hold [37, 65, 66, 42, 43]. 
Elliptic models. At the surface of maximal expansion {Hq =0), we have from (46) 

f= 1 + 3(A[,''^ - A^'^^). (72) 

In the limit a — )■ with t — )■ ibb we obtain the same expression (69) as in the 
hyperbolic case in the same limit, but in the limit a — >■ with t — )• teon we have 
'Hq — oo, hence we obtain 

f «l + 3(A('')-Af)) + [^^, (73) 

where 

l-r£'/£~^ V ° 2 ^0 ) + 1 - r£'/£ ' 

ft > implies M' — ME' /EE > 0, this means that unlike in the LTB model we cannot have M rj 
constant. 
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follows from (57c). Hence, considering (69) and (72)-(74), we have the following 
conditions for F > in terms of Szekeres initial fluctuations 

1 + 3(A^'') - A^'^)) > 0, A^'') - ^A^'^) > 0, necessary, (75a) 

A^,''^ > -1, ^ < 0, ^ > 0, necessary and suflicient, (756) 

1 — rt' /c 1 — rt' /c 

As in the paraboUc and hyperbolic cases, these conditions reduce to the Hellaby- 

Lake conditions in the LTB limit £' = 0, A^^' = [37, 42, 43, 65, 66]. 

9. Comparison with LTB models. 

Since LTB models arc well known inhomogeneous cosmological models that have been 
frequently utilized, it is useful and desirable to compare them with the quasi-spherical 
Szekeres models, which provide a straightforward non-spherical generalization. 

The evolution equations (31a)-(31(i) are identical to those of an LTB model [37], 
save for the fact that A'^''^ and A^^' depend on {x,y), though this dependence only 
needs to be prescribed as part of the initial conditions through the term 1 — r£' /£, 
with £ given by (5) or £' /£ by (43) in spherical coordinates. Likewise, the forms of 
the QL scalars Pq, fCg, Hg and fig given by the scaling laws (48a)-(50) are identical to 
those of an LTB model described by the q-scalars and their fluctuations [37, 42, 43] 
(notice that the determination of the scale factor a in (53)-(55) does not require that 
£ is specified). 

Since the q-scalars are common to both LTB and Szekeres models, the difference 
between these models only enters in the fluctuations through the term 1 — r£' /£, as 
can be appreciated from the relation between initial fluctuations Aq'^'' and 6q^^ of 
Szekeres and LTB models given in (39a)-(396): 



Considering the fact that any LTB model can be characterized as a unique solution 
of the system (31a) (Sid) for initial value functions {pgo, ICqo, cJq'''', '^o'^^}, then each 
LTB model can be associated with a Szekeres model by the following transformation 
in the space of initial conditions: 

4^)^A^^\ so that 5(^)^->A(^), (77) 

with A(^) for A = m,k,n given by (486), (48c), (51) and (35). The transformation 
(77) simply requires modifying the initial fluctuations of an arbitrary LTB model by 
choosing (besides pgo, JCqo) the three extra free functions {S*, P, Q} to construct the 
term £' /£ in (43). Szekeres models obtained in this manner form a 3-parameter 
class of models associated with a unique LTB model that follows from the solution 
of the same evolution equations [i.e. (31a)-(31(i)) but with the modified initial 
conditions {pgo, /Cgo, Aq''"', Aq'^''}. Conversely, any quasi-spherical Szekeres model 
can be mapped to a unique LTB model with the same QL scalars pg, Kg, Tig and fig 
(given by (48a)-(50)) and fiuctuations transformed by (77). 

In particular, the relation between a given LTB model and the 3-parameter class 
of associated Szekeres models can be understood in terms of a perturbative approach 
if we choose the free functions {S, P, Q} so that r£'/£ 4C 1, which implies: 

A^"^-fl + ^)4^ (78) 
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Figure 1. Evolution of F for a hyperbolic model with negligible dipole (r£' /£ !V 
0), and with ilg = 0.3 and -Hq = 70 km s"! Mpc'^ (see (34)) Upper left 
(a) S^'^^ = -0.66, upper right (b) 5^'^' = -0.1, lower left (c) ^'^^ = 0.5, 
lower right (d) (5q = 1. In each panel the curves from the top to bottom: 
5^''' = 1, 0.5, 0, -0.5, and -0.99. 



SO that the initial fluctuations Aq ' take the form of perturbations of the LTB 
fluctuations. Under these conditions, we can examine Szekeres models that are almost 
LTB, with perturbative deviations from spherically symmetric. It is important to 
remark that quasi plane and quasi-hyperbolic Szekeres models relate to dust solutions 
with plane and pseudo-spherical symmetry in the same manner as quasi-spherical 
models relate to spherically symmetric LTB solutions (which we discussed in this 
section). 

10. Numerical example: Growth of the dipole distribution. 

In order to illustrate how the theoretical framework that we have presented works in 
practice, we examine the dipole evolution that marks the deviation of a quasi-spherical 
Szekeres model from spherical symmetry. In particular, we address the question of the 
stability of spherical symmetry with respect to dipole perturbations. Let us consider 
a hyperbolic {K < 0) model with initially small deviation from spherical symmetry, 
i.e. with r£' /£ w (see (78)). The evolution of this model for a comoving surface of 
fixed r is calculated as follows: 

(i) First we choose an FLRW background that can be identified in the asymptotic 
limit r — )• oo of f2q and "Hq in the slice t = to. Using (34), we choose as background 
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Figure 2. Evolution of F (from to till tcoii) for a elliptic model with negligible 
dipole (rS'/S ^ 0), and with Ug = 0.3 and Hq = 70 km s"! Mpc"! (see (34)) 
left (a) <5g'^' = —0.66, right (b) 6^^^ = —0.1. In each panel the curves from the 
top to bottom: 4''' = 1,0.5,0, -0.5, and -0.99. 



quantities fio = ^qoo = 0.3 and Hq = Hgoo = 70 km s ^ Mpc ^. f 

(ii) We choose the initial value functions pqo and K.qo at a fixed r with t = to 
corresponding to the present cosmic time (computed from (54) considering present 
day values for Hq and flo given by the background values mentioned in point (i) 
above). The initial monopole perturbations Sjf^ and Sq^^ are computed from 
these initial condition from (76). 

(iii) We consider a small and negligible dipole perturbation. One of such choices is: 

P' = = Q', 5 = r" with a ^ 0.001. 

(iv) Knowing (5q''\ Sq'^'' and £'/£ we find Ag''-' and Aq'^-' from from (76). 

(v) As seen from (486), under these conditions the density perturbation evolves only 
through the time evolution of F, which is given by and (58a) or (586). Notice 
that the form of A^''' is not only determined by F, but also by the presence of 
r£' /£ in Aq''\ which, however, does not depend on time. 

Since initially Aq''^ ft! ^q''^ (i.e. r£'/£ « 0) thus the only possibility for a 
large departure from the spherical symmetry (large dipole-like fluctuations) is when 
r w r£'/£. 

The evolution of F is presented in Fig. 1. As seen in the cases that were 
considered, F tends to an asymptotic value Too- Thus, we recover the know fact from 
the evolution of ever-expanding hyperbolic LTB models (i.e. K < 0), that density 
perturbations freeze in the time asymptotic range [16], but only if F^o « we can 
have an asymptotically large dipole variation. If Too » r£' /£ then the deviation from 
spherical symmetry is negligible. As can be seen from Fig. 1, only when (5q < we 
have Foo < 1- Thus, if the initial deviation from spherical symmetry is small then the 
spherical shape is stable as long as we exclude models that are close to shell crossing 
singularities (i.e. Sq'^^ « —2/3). Wc; can also obtain this result by looking at the full 
analytic expression for the density perturbation A'''^ from the scaling law (486) and 

t Since the examples we are considering are meant for illustrative purposes, these values simply allow 
us to specify a model and do not correspond to an actual model of some kind of a realistic structure. 
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the form of f in (586) (and eliminating A{,''\ A^'^' in terms of 6lf \ ^'^^ with (76)): 

3 (4") - |<5f [n,it t,,) - f ] + rn,t[, 
A^P^ = - ^— - (79) 

1 - r£'/£ + 3(4") - 4^^)) - 3 (4") - f 5f )) H,{t t,,) ~ rH,C 

In the Hmit i — >■ oo (which is equivalent to a — >■ cxd) we obtain 'Hg — ^ from (49) and 
'^iq{t — ibb) 1 from (54) and (70). Hence, we obtain in this limit 

1 - rS'/£ + §4'=) V « y 

which shows the value at which density perturbations freeze asymptotically, and 
also the fact that the asymptotic effects of the deviation from spherical symmetry 
(contained in the term r£' /£) remain small if chosen to be small at t = to. 

We now consider another similar numerical example. This time, however, we 
take an elliptic background, thus we use the analytic solution (55) and consider a 
value rtqo = 1.5 > 1 for a fixed r. The results are presented in Fig. 2. The evolution 
is calculated from the initial instant to (present cosmic time) all the way to the big 
crunch t = tcou- As can be seen in both panncls of the figure, there arc examples in 
which a shell crossing singularity occurs (F = 0) , and it is easy to verify that for these 
values of Sq'^ and Sq^^ the regularity conditions (75a)-(756) are not satisfied. If we 
demand that the models are free of shell crossings, then F ^ must not occur, hence 
F must increase monotonically as dust layers expand for early times and collapse 
for late times (there is not symmetry with respect to t = tbou„ce where Hq = 0). 
Thus, the density perturbation A^''^ monotonically decreases, and as F ^ oo, we have 
AW -1. Hence, a generic perturbation within the elliptic background will either 
decrease (A^''^ — >■ —1) or increase and eventually lead to a shell crossing. 

The effects of the deviation from spherical symmetry also remain small (if 
originally small at t = to)- This can be seen looking at the analytic form of A^''^ 
in (79) in the limit t t^ou- Since 'Hg —oo as t toon, then an asymptotic 
expansion of A^^'^ around — \Hq\ yields: 

- -1 + '^if^i^/f^^ + o{\nqn, (81) 

where wc used (74). Evidently, the effects of the deviation from spherical symmetry 
(the term r£' /£) enter as a correction of order and thus remain small if chosen 

initially to be small. The limit A^'') ^ — 1 as t — > ieoii rnay seem strange, as one would 
associate a diverging density contrast near the collapsing singularity. This limiting 
value follows from (486) and from the fact that F and F diverge as Hq — > — oo (which 
occurs as t ^ icon)- Notice that the A^"^' arc NOT "contrast" perturbations, but 
have a more complicated interpretation related to the time evolution of the gradients 
A' and A'^ through their definition in (24). Therefore, there is no reason for A('') to 
diverge at t^^u- The limit A^''^ — )• —1 could imply Schwarzschild vacuum conditions if 
p ^ with pq > [42, 41], but in the collapsing regime it reflects the fact that the 
ratio p/pq vanishes as i — >■ t^ou, but with both densities diverging in this limit. 

The following conclusion arises from the discussion above: if we consider only 
initial perturbations and demand absence of shell crossings together with an almost 
spherical shape with negligible dipole-like departure from spherical symmetry, then 
the spherical shape is conserved. Bearing in mind this conclusion, the following 
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question arises: If spherical symmetry is a stable property, then why do non- 
symmetrical structmes arc present in the Universe? The reason is that even if cosmic 
structures (voids, clusters, superclusters) evolved from small fluctuations that were 
present at the last scattering time, these small fluctuations did not have to be almost 
spherical, nor the conditions of avoidance of shell crossings had to be satisfied initially. 
In fact, the only case in which an initially large deviation from spherical symmetry 
allows us to permanently (not temporally - see insets in Fig. 1) dissolve the dipole 
occurs only when the decaying modes are present. One example of this is a parabolic 
model with zero curvature perturbations, i.e. /Cgo = = Sq^\ Then as seen from 
(58a) foo = 1 + A^''^ and so A^'') 0. 

11. Final discussion and further work. 

We have introduced for Szekeres models of class I a set of new coordinate independent 

representations of scalar variables consisting of the q-scalars and their fluctuations. 
We have shown throughout the article that these variables completely determine 
the dynamics of the models, either in terms of intuitively appealing analytically 
expressions (section 7) or through fluid flow evolution equations (section 5) that 
can be handled as ordinary differential equations (ODE's) subjected to algebraic 
constraints. We have also shown that by applying these constraints we can construct 
various equivalent representations of the new variables, some of which can be related 
to observational parameters such as the Hubble factor H and Q (sections 5.2 and 7.3). 
Various related theoretical issues have been discussed in detail: 

— the relation between the q-scalars and averages of covariant scalars (section 4), 

— initial conditions for the fluid flow evolution equations (section 6), 

— regularity conditions to avoid shell crossings (section 8), 

— as well as a comparison with LTB models (section 9), 

— we have also used the new variables (section 10) to examine the preservation of 
nearly spherical initial conditions in late time regimes, and showed that spherical 
symmetry is stable against small dipole-like perturbations. 

We list below the main advantage of these variables over the traditional ones, as 
well as potential applications to be undertaken in future articles: 

Numerical work. The fluid flow evolutionary equations that were obtained in 
section 5 (in any given representation) form a system of four PDE's that can 
be handled effectively as ordinary differential equations (ODE's). One possible 
representation is that given by (31a)-(31d), consisting of the q-scalars associated 
with the density and expansion scalars and their fluctuations with respect to the 
local scalars. In another representation (in (36a)-(36rf)) the four variables are the 
q-scalars associated with the Hubble and Omega factors and their fluctuations. 
In either representation, the system depends on five free parameters that convey 
the effects of spherical and non-spherical inhomogeneity and are specified as 
initial conditions. The spacelike constraints for these system are not PDE's, 
but algebraic equations. This represents an important advantage over the fluid 
flow evolution equations in terms of the local covariant scalars, which need to 
be handled as PDE's because the constraints are PDE's that couple in a non- 
trivial way to the time derivatives (as in the fiuid flow systems in [53, 54]). 
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These evolution equations provide a nice and elegant (and simplified) approach to 
numerical work with Szckcrcs models, with an enormous potential for applications 

cither in theoretical studies or in fitting observations. 

Theoretical work. By employing an initial value formalism based on the q-scalars 
and their fluctuations, we can 

• study the propagation of any given set of initial conditions, either analytically 
or numerically. We can compare this propagation with that of initial 
conditions in spherical LTB models. In particular, we can examine the 
stability of initial spherical shapes (LTB model) against small (dipole like) 
perturbations - see section 10. 

• extend previous theoretical results obtained by means of these variables in 
LTB solutions to the Szckcrcs models. In particular, we aim in future articles 
to generalize previous work dealing with important features of the LTB 
models: averaging inhomogeneities [36, 38, 39, 40, 41], radial asymptotics 
[42] , evolution of radial profiles [43] , and dynamical systems analysis [44, 45] . 

• extend the work done on the connection with non-linear perturbations on 
a FLRW background [35] that compares the fluctuations between local and 
QL variables with exact perturbations in a FLRW background, for example 
identifying and studying the evolution of exact quantities that generalize the 
growing and decaying modes of the theory of linear perturbation of dust 
sources. 

Wc consider that the new variables that we have introduced not only provide a 
deeper theoretical understanding of the Szekeres models, but is (at the same time) 
more intuitive than the study of the models in the traditional variables. This formalism 
has an enormous potential for exploring the effects of non-spherical inhomogeneity and 
non-linear perturbations, it may also allow for a more efficient utilization of Szekeres 
solutions in for the study of cosmic inhomogeneities (including void models) to test 
cosmological observations. We are currently undertaking further work on these lines 
that we expect to submit in the near future. 

Appendix A. Averaging 

As mentioned in Section 3, the integral definition of the q scalars is equivalent to a 
proper volume averaged integral with weight factor T. The standard proper volume 
averaging (weight factor = 1) for the quasi-spherical Szekeres models was studied 
in [56]. Wc show in this appendix that, save for some qualitative differences, the 
results and the approach of [56] remain valid for the quasi-hyperbolic and quasi- 
planar models, leading to the expressions (19), (20) and (21). 

First, we remark that the area of a 2-surface of constant t and r in the quasi- 
hyperbolic and quasi planar models may be infinite. However, this is not problematic, 
as the surface of the domain Su cancels out. Second, a location r — that can 
be identified as an "origin" only exists for quasi-spherical models - in the quasi- 
hyperbolic model r cannot be equal to zero, and in the quasiplane r can only 
asymptotically approach the origin, r — >■ [15]. However, one can always consider a 
domain V is centered around r = even if this point does not belong to the manifold. 

For the purpose of averaging it is more convenient to adopt a pair of complex 
conjugate coordinates 

C = x + iy, C = x-iy, (A.l) 
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(A.2) 



The 3 vohimc associated with the quasi-local average of the domain V centered at 
the origin is (the weight factor = y/e — K, see also (18)) 



jdrjj dCdC £Y'Y^ = j'^^'j j ~ h 



= / dr 



dCdC , 1 „3 d 
£2 ^ 2 dr 



£2 y 



(A.3) 



where Tc is the lower limit of integration (only for the quasi-spherical and quasi-planar 
model Tc = 0), dCdC/£2 is the metric of a unit sphere/plane/hyperboloid and 

dCdC 



£2 



which does not depend on r (for the quasi-spherical model St> = 47r). Thus 

Vq = S-D [ drR^R' = S-dT-d (A.4) 



Note that R{rc) can be equal zero [15], so Tp = (l/3)i?|,, even if for the quasi- 
hyperbolic and quasi-planar model r cannot be equal to 0. The q-density pq is 



Pg = - I dr I I dCdC SY'Y'p 



£'\ 2M'-QM£'/£ 
J) i?2 {R' - R£'/£) 



or 1 
^ I drM' + 



. drM — 
35x,Tp ./ dr 



dcdc\ 

£2 J 



1 /drM' = ^ 
TvJ Rl 



The q-scalar Og dual to the local expansion © is 



(A.5) 



r-D 



e, = -idr 



dCdC o2 /^p/ _ i?' + 2RR'/R - 3R£'/£ 

R' - R£'/£ 



£2 '-R\R'-R^ 



, „9 „/ / R' ^ R 



Rt> 
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Finally, the q-scalar ^TZq dual to ^TZ is 



r-D 



= / dr{RKy = 6^. (A.7) 

ra 

Note that (A.5), (A.6) and (A.7) exactly coincide with (19), (20) and (21) at the 
domain boundary marked by r. 

Appendix B. CovariEint meaning of the q— scalars and their fluctuations. 

It is straightforward to show that the q-scalars pq, Hq, ICq, flq and their fluctuations 
are coordinate indcipcndciiit quantities expressible in terms of curvature invariants. 
The q-scalar pg and the fluctuation A^''^ take the form 

8npq = 6^,-11, A(^) = ^^^, (B.l) 

where TZ = Sirp is the 4 -dimensional Ricci scalar and ^['2 = ^VV is the only nonzero 
Weyl curvature invariant in a Newman-Penrose tetrad. The fact that the ratio of 
Weyl to Riemann curvature scalars can provide a measure of inhomogeneity by means 
of suitably defined density fluctuations in LTB models has been already highlighted 
in [67]. Expressions similar to (B.l) follow from (29&) for T-Lq and A*^^); 

nq = n + i:, a(«) = -3^^, (b.2) 

where H = 38 = h^Vbu"' is the expansion scalar and S is the scalar in (11) associated 
with the shear tensor. While S is already given in a coordinate independent form 
in (11), it is easier to provide an invariant characterization of it as the independent 
eigenvalue of the shear tensor fi" 5. Notice that the form of A^^^ above is directly 
related to the scalar ratio S/0 or the quadratic form aab<7°'^/'d^ = 65]^/0^ which 
provide coordinate independent measures of the ratio of anisotropic vs isotropic 
local expansion velocities. Coordinate independent forms for the remaining q-scalars 
ICq, rig and fluctuations A^'^^ A^^^^ follow by applying the constraints (32)-(33) and 
(34)-(35) to (B.l) and (B.2). 

Appendix C. Functions of q— scalars are q— scalars. 

Let i/j = tjj{Aq,Bq) be a smooth function of Aq and Bq complying with (18), then 
considering (22) and (24), ip = (j)q is the q-scalar for which we identify A'*") and </> 
by: 

Ip dAq Xjj ABq 



This property was used to define A*^^^ and f2 in (35), since — ^q(^pq^T~iq^ in 
(34). However, as opposed to the local scalars p, H, fC, the physical and geometric 
interpretation of the "local" scalar SI is not clear. 
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Appendix D. Models with "closed" and "wormhole" topologies. 

The coordinate choice (37) cannot be used if the shces ^T[t] are "closed" 
(homeomorphic to or are "wormholes" (homeomorphic to either x M or §^ x S^). 
In the former case there arc two worldhnes where Rq = R = 0, which generahze the 
closed spherical models with two symmetry centers and in the latter Rq has no zeroes. 
In both cases Rq must have a zero. Hence, we choose: 

Ro = iof{r), (D.l) 

where io is an arbitrary length scale and /(r) is a dimensionless function with the 
appropriate properties. For closed models / can be a sine-type of function, while for 
the wormhole case it can be either cosh or sec. The coordinate choice (D.l implies 
the following replacement 

5V5_dln^ 

S ^ /'//"din/ ^"-'^ 
hence, /' and £' (and thus S', P' , Q') must have same order common zeroes. 
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